% This practical will work with a single subject's data from an emotional
% faces experiment (Elekta Neuromag data). You can get the data from:
% https://sites.google.com/site/ohbaosl/practicals/practical-data/emotional-face-processing-elekta-neuromag-data
% 
% Work your way through the script cell by cell using the supplied dataset.
% As well as following the instructions below, make sure that you read all
% of the comments (indicated by %), as these explain what each step is
% doing. Note that you can run a cell (marked by %%) using the ?Cell? drop
% down menu on the Matlab GUI.    

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
% set this to where you have downloaded OSL and the practical data:
practical_dir='/home/mwoolrich/Desktop'; 
osldir=[practical_dir '/osl1.3.1'];    

practical_dir='/Users/woolrich';
osldir=[practical_dir '/homedir/matlab/osl1.3.1'];    

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

% directory where the data is:
%workingdir=[practical_dir '/faces_subject1_data']; % directory where the data is
workingdir=[practical_dir '/homedir/matlab/osl_testdata_dir/faces_subject1_data']; % directory where the data is

cmd = ['mkdir ' workingdir]; if ~exist(workingdir, 'dir'), unix(cmd); end % make dir to put the results in

clear spm_files_continuous spm_files_epoched structural_files;

% set up a list of SPM MEEG object file names (we only have one here)
spm_files_continuous{1}=[workingdir '/spm8_meg1.mat'];
spm_files_epoched{1}=[workingdir '/espm8_meg1.mat'];

% set up a list of mris
structural_files{1}=[workingdir '/structurals/struct1.nii'];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SETUP BEAMFORMER AND FIRST-LEVEL GLM OAT
% In this section we will do a wholebrain beamformer, followed by a trial-wise
% GLM that will correspond to a comparison of the ERFs for the different
% conditions.

%% SETUP THE OAT:
oat=[];
%oat.source_recon.D_continuous=spm_files_continuous;
oat.source_recon.D_epoched=spm_files_epoched;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
do_tf=1;
if do_tf,
    oat.source_recon.freq_range=[1 30]; % frequency range in Hz
    oat.source_recon.time_range=[-0.2 0.4];
    
    oat.first_level.time_range=[-0.15 0.35];
    oat.first_level.tf_method='hilbert';
    oat.first_level.tf_hilbert_freq_res=19;
    
    oat.first_level.tf_hilbert_do_bandpass_for_single_freq=0;
    %oat.first_level.tf_hilbert_freq_ranges=[4 8; 8 13; 13 30; 1 48];
    oat.first_level.tf_hilbert_freq_ranges=[1 48];
    oat.first_level.tf_num_freqs = size(oat.source_recon.freq_range,1);
    
    oat.first_level.tf_time_step = 0.015;
else
    oat.source_recon.freq_range=[1 48]; % frequency range in Hz
    oat.source_recon.time_range=[-0.2 0.4];
end;
%oat.source_recon.time_range=[-0.1 0.3];
oat.source_recon.method='beamform';
%oat.source_recon.method='mne';

oat.source_recon.gridstep=8; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.mri=structural_files;
oat.source_recon.modalities={'MEGPLANAR','MEGMAG'};
%oat.source_recon.modalities={'MEGPLANAR'};
oat.source_recon.bandstop_filter_mains=0;

%oat.source_recon.forward_meg='MEG Local Spheres';
oat.source_recon.forward_meg='Single Shell';

oat.source_recon.work_in_pca_subspace=0;
oat.source_recon.force_pca_dim=1;
oat.source_recon.pca_dim=55;

do_hmm=1;
if(do_hmm)
    oat.source_recon.hmm_num_states=3;
    oat.source_recon.hmm_num_starts=1;
    oat.source_recon.hmm_pca_dim=30;
    oat.first_level.hmm_do_glm_statewise=0;
else
    oat.first_level.hmm_do_glm_statewise=0;
end;

% design_matrix_summary is a parsimonious description of the design matrix.
% It contains values design_matrix_summary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to create the (num_regressors x num_trials) design matrix:
design_matrix_summary={};
design_matrix_summary{1}=[1 0 0 0];design_matrix_summary{2}=[0 1 0 0];design_matrix_summary{3}=[0 0 1 0];design_matrix_summary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=design_matrix_summary;

% Alternatively you can specify trial-wise subject-specific design matrices
% as ASCII files
if(0),
    % Need to create design matrix without rejects. D.pickconditions auto
    % excludes rejected trials:
    De=spm_eeg_load(spm_files_epoched{1});
    Ntrials=De.ntrials;
    X=zeros(Ntrials,length(oat.source_recon.conditions));
    tot=0;
    for ii=1:length(oat.source_recon.conditions),
        num=length(De.pickconditions(oat.source_recon.conditions{ii}));
        X(tot+1:tot+num,ii)=1;
        tot=tot+num;
    end;
    X=X(1:tot,:);
    
    % save as ASCII file:
    save([workingdir '/subject1_desmat.txt'],'X','-ascii');
    design_matrix_summary={};
    design_matrix_summary{1}=[workingdir '/subject1_desmat.txt'];
    
    try,rmfield(oat.first_level,'trial_rejects');catch, end;
    oat.first_level.design_matrix_summary=design_matrix_summary;
end;

tname=[regexprep(oat.source_recon.forward_meg, ' ', '') '_hmm' num2str(do_hmm) '_' num2str(oat.source_recon.method) '_gsw' num2str(oat.first_level.hmm_do_glm_statewise) ];
oat.source_recon.dirname=[spm_files_continuous{1} '_wideband_' tname '_' osl_version];

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast_name={};
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';
oat.first_level.cope_type='acope';

oat = osl_check_oat(oat);

%% RUN THE OAT:

oat.to_do=[1 1 0 0];
oat.first_level.name=['wholebrain_first_level'];
oat = osl_run_oat(oat);

res=osl_load_oat_results(oat,oat.source_recon.results_fnames{1});

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% view the GLM design matrix (NOTE that column 1 is motorbikes, columns 2-4 are faces)
figure;imagesc(stats.x);title('GLM design matrix');xlabel('regressor no.');ylabel('trial no.');

% results = osl_get_recon_timecourses( osl_load_oat_results(oat,oat.source_recon.results_fnames{1});, mask_fname )
% osl_oat_plot_hmm_states( oat )

%% OUTPUT SUBJECT'S NIFTII FILES
% Having run the GLM on our source space data, we would like to inspect the
% results for our single subject. 
% We can do this by saving the contrast of parameter estimates (COPEs) and 
% t-statistics for each of our contrasts to NIFTI images.

S2=[];
S2.oat=oat;
%S2.stats=stats;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[3]; % list of first level contrasts to output
S2.resamp_gridstep=oat.source_recon.gridstep;
[statsdir,times,count]=osl_save_nii_stats(S2);

% VIEW NIFTII RESULTS IN FSLVIEW
% We can now view the nifti images containing our GLM results in FSL, here
% we are runnin fslview from the matlab command line, but you do not need
% to - you can run it from the UNIX command line instead.

mni_brain=[OSLDIR '/std_masks/MNI152_T1_' num2str(S2.resamp_gridstep) 'mm_brain']; 

% INSPECT THE RESULTS OF A CONTRAST IN FSLVIEW. Recall:
%S2.contrast{1}=[3 0 0 0]'; % motorbikes
%S2.contrast{2}=[0 1 1 1]'; % faces
%S2.contrast{3}=[-3 1 1 1]'; % faces-motorbikes

con=3;
runcmd(['fslview ' mni_brain ' ' [statsdir '/cope' num2str(con) '_' num2str(S2.resamp_gridstep) 'mm'] ' ' [statsdir '/pseudo_zstat_var_' num2str(S2.resamp_gridstep) 'mm'] ' ' [statsdir '/tstat' num2str(con) '_' num2str(S2.resamp_gridstep) 'mm']  ' ' [statsdir '/tstat' num2str(con) '_mip_' num2str(S2.resamp_gridstep) 'mm'] ' &']);

